#include "tools.h"
#include <iostream>

#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;} 
#define TINY 1.0e-20 
#define MAXITS 1000 
#define EPS 1.0e-5
#define NR_END 1



void skip_line( FILE * fp ) { // Skip a text line in an input file
  char ch;
  do ch=getc(fp);
  while(ch!='\n');	
}
	
double next( FILE * fp ) { // Reads a double value at begining of a line of a text file
  double buffer;
  skip_line(fp);
  fscanf_s(fp,"%lf",&buffer);
  return buffer;
}

void nrerror( char * str) { // Print a error message and exit the program.
  printf("Error: %s\n",str);
  printf("Program exits\n");
  exit(1);
}


void rest() { // Wait for you to press the enter key on the keyboard 
  char ch='x';
  printf("Return to continue...or q to quit  ");
  while(ch!='\n') {
      if(ch=='q') exit(1);
      else scanf_s("%c",&ch);
  }
  
}

void noMoreMemory() {  // Function that handle the failure of the new() command.
  
  printf("Unfortunately, the dynamic allocation of memory has failed!\n");
  printf("This means you run out of memory and your program is going to stop.\n");
  printf("Now, it could be a bug(I have forgotten to release memory) or a system shortage\n");
  printf("I can't help you with that\n");
  printf("ABORT!\n");
  abort();
  
}



/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

double *vector(int nl, int nh) {

/* allocate a double vector with subscript range v[nl..nh] */

	double *v;
	
	v = new double [nh-nl+1+NR_END]; 
	if (!v) nrerror("allocation failure in vector()"); 
	return v-nl+NR_END; 	
} 

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
int *ivector(int nl, int nh) { 

/* allocate an int vector with subscript range v[nl..nh] */ 

	int *v; 

	v  = new int [nh-nl+1+NR_END]; 
	if (!v) nrerror("allocation failure in ivector()"); 
	return v-nl+NR_END;
}


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
double **matrix(int nrl, int nrh, int ncl, int nch) {

/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */

	int i, nrow,ncol; 
	double **m; /* allocate pointers to rows */

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	m = new double* [nrow+NR_END];
	if (!m) nrerror("allocation failure 1 in matrix()"); 
	m += NR_END; 
	m -= nrl; /* allocate rows and set pointers to them */ 
	m[nrl] = new double [nrow*ncol+NR_END];
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); 
	m[nrl] += NR_END; 
	m[nrl] -= ncl; 
	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ 
	return m; 

} 

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
double ***tensor(int nrl, int nrh, int ncl, int nch, int ndl, int ndh) {

/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */

	int i,j,nrow,ncol,ndep; 
	double ***t; 

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	ndep=ndh-ndl+1;
	t = new double** [nrow+NR_END]; 
	if (!t) nrerror("allocation failure 1 in tensor()"); 
	t += NR_END; 
	t -= nrl; /* allocate pointers to rows and set pointers to them */ 
	t[nrl] = new double * [nrow*ncol+NR_END]; 
	if (!t[nrl]) nrerror("allocation failure 2 in tensor()"); 
	t[nrl] += NR_END; 
	t[nrl] -= ncl; /* allocate rows and set pointers to them */ 
	t[nrl][ncl] = new double [nrow*ncol*ndep+NR_END]; 
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in tensor()");
	t[nrl][ncl] += NR_END; 
	t[nrl][ncl] -= ndl; 
	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep; 
	for(i=nrl+1;i<=nrh;i++) { 
		t[i]=t[i-1]+ncol; 
		t[i][ncl]=t[i-1][ncl]+ncol*ndep; 
		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep; 
	} /* return pointer to array of pointers to rows */ 
	return t; 

}


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* free a double vector allocated with vector() */

void free_vector(double *v, int nl, int nr) {delete [] (v+nl-NR_END);}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* free an int vector allocated with ivector() */

void free_ivector(int *v, int nl, int nr) { delete [] (v+nl-NR_END); }


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
void free_matrix(double **m, int nrl, int nrh, int ncl, int nch) {

/* free a double matrix allocated by matrix() */

	delete [] (m[nrl]+ncl-NR_END); 
	delete [] (m+nrl-NR_END); 
} 

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void free_tensor(double ***t, int nrl, int nrh, int ncl, int nch, int ndl, int ndh) {

/* free a double f3tensor allocated by f3tensor() */
	
	delete [] (t[nrl][ncl]+ndl-NR_END); 
	delete [] (t[nrl]+ncl-NR_END); 
	delete [] (t+nrl-NR_END); 

}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void bandec(double **a, int n, int m1, int m2, double **al, int indx[], double *d) {

  /*Given an nXn band diagonal matrix A with m1 subdiagonal rows and m2 
    superdiagonal rows, compactly stored in the array a[1..n][1..m1+m2+1] 
    this routine constructs an LU decomposition of a rowwise permutation of A.
    The upper triangular matrix replaces a, while the lower triangular matrix 
    is returned in al[1..n][1..m1]. indx[1..n] is an output vector which records 
    the row permutation effected by the partial pivoting; d is output as   1 
    depending on whether the number of row interchanges was even or odd, respectively.
    This routine is used in combination with banbks to solve band-diagonal sets of equations.*/

  int i,j,k,l;
  int mm;
  double dum;

  mm=m1+m2+1;
  l=m1;
  for (i=1;i<=m1;i++) { /*Rearrange the storage a bit.*/
    for (j=m1+2-i;j<=mm;j++) a[i][j-l]=a[i][j]; 
    l--; 
    for (j=mm-l;j<=mm;j++) a[i][j]=0.0; 
  } 
  *d=1.0; 
  l=m1; 
  for (k=1;k<=n;k++) { /*For each row... */
    dum=a[k][1]; 
    i=k; 
    if (l < n) l++; 
    for (j=k+1;j<=l;j++) {	/*Find the pivot element.*/
      if (fabs(a[j][1]) > fabs(dum)) { 
	dum=a[j][1]; 
	i=j; 
      } 
    } 
    indx[k]=i;
    if (dum == 0.0) a[k][1]=TINY; 
    /*Matrix is algorithmically singular, but proceed anyway with TINY pivot 
      (desirable in some applications). */
    if (i != k) { /*Interchange rows.*/
      *d = -(*d); 
      for (j=1;j<=mm;j++) SWAP(a[k][j],a[i][j]) 
			    } 
    for (i=k+1;i<=l;i++) { /*Do the elimination.*/
      dum=a[i][1]/a[k][1];
      al[k][i-k]=dum; 
      for (j=2;j<=mm;j++) a[i][j-1]=a[i][j]-dum*a[k][j]; 
      a[i][mm]=0.0; 
    } 
  } 

}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */


void banbks(double **a, int n, int m1, int m2, double **al, int indx[], double b[]) {

/*Given the arrays a, al, and indx as returned from bandec, and given a right-hand 
side vector b[1..n], solves the band diagonal linear equations A   x = b. 
The solution vector x overwrites b[1..n]. The other input arrays are not modified, 
and can be left in place for successive calls with different right-hand sides. */

	int i,k,l;
	int mm;
	double dum;
	
	mm=m1+m2+1; 
	l=m1; 
	for (k=1;k<=n;k++) { /*Forward substitution, unscrambling the permuted rows as we go.*/
		i=indx[k]; 
		if (i != k) SWAP(b[k],b[i]); 
		if (l < n) l++;
		for (i=k+1;i<=l;i++) b[i] -= al[k][i-k]*b[k]; 
	} 
	l=1;
	for (i=n;i>=1;i--) { /*Backsubstitution.*/
		dum=b[i]; 
		for (k=2;k<=l;k++) dum -= a[i][k]*b[k+i-1];
		b[i]=dum/a[i][1];
		if (l < mm) l++;
	} 

}

/* ------------------------------------------------------------ */

void tridag(double a[], double b[], double c[], double r[], double u[], unsigned long n) {

  /*Solves for a vector u[1..n] the tridiagonal linear set given by equation (2.4.1). 
    a[1..n], b[1..n], c[1..n], and r[1..n] are input vectors and are not modified. */

  unsigned long j;
  double bet,*gam;

  gam=vector(1,n); /*One vector of workspace, gam is needed. */
  if (b[1] == 0.0) nrerror("Error 1 in tridag"); 
  /*If this happens then you should rewrite your equations as a set of order N-1 with u2
    trivially eliminated*/
  u[1]=r[1]/(bet=b[1]);
  for (j=2;j<=n;j++) { /*Decomposition and forward substitution.*/
    gam[j]=c[j-1]/bet; 
    bet=b[j]-a[j]*gam[j];
    if (bet == 0.0) nrerror("Error 2 in tridag"); /*Algorithm fails; see below.*/
    u[j]=(r[j]-a[j]*u[j-1])/bet;
  } 
  for (j=(n-1);j>=1;j--) u[j] -= gam[j+1]*u[j+1]; /*Backsubstitution.*/
  free_vector(gam,1,n); 
}

/* ------------------------------------------------------------ */

void lubksb(double **a, int n, int *indx, double b[]) {

  /*Solves the set of n linear equations A   X = B. Here a[1..n][1..n] is input,
    not as the matrix A but rather as its LU decomposition, 
    determined by the routine ludcmp. indx[1..n] is input as the permutation 
    vector returned by ludcmp. b[1..n] is input as the right-hand side vector B, 
    and returns with the solution vector X. a, n, and indx are not modified by this 
    routine and can be left in place for successive calls with different right-hand sides b.
    This routine takes into account the possibility that b will begin with many zero
    elements, so it is efficient for use in matrix inversion. */
  
  int i,ii=0,ip,j;
  double sum;
  
  for (i=1;i<=n;i++) { 
    /*When ii is set to a positive value, it will become the index of the first 
      nonvanishing element of b. We now do the forward substitution. The only new 
      wrinkle is to unscramble the permutation as we go. */
  ip=indx[i]; 
  sum=b[ip]; 
  b[ip]=b[i]; 
  if (ii) 
    for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; 
  else if (sum) ii=i; 
  /*A nonzero element was encountered, so from now on we will have to do the sums in the loop above.*/ 
  b[i]=sum; 
  } 
  for (i=n;i>=1;i--) { /*Now we do the backsubstitution.*/
    sum=b[i]; 
    for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; 
    b[i]=sum/a[i][i]; /*Store a component of the solution vector X. */
  } 
}

/* ------------------------------------------------------------ */

void ludcmp(double **a, int n, int *indx, double *d)  {

/*Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition
 of a rowwise permutation of itself. a and n are input. a is output, arranged as in
 equation (2.3.14) above; indx[1..n] is an output vector that records the row
 permutation e ected by the partial pivoting; d is output as   1 depending on
 whether the number of row interchanges was even or odd, respectively. 
 This routine is used in combination with lubksb to solve linear equations or
 invert a matrix. */

  int i,imax,j,k;
  double big,dum,sum,temp;
  double *vv; /*vv stores the implicit scaling of each row.*/

  vv=vector(1,n);
  *d=1.0; /*No row interchanges yet.*/
  for (i=1;i<=n;i++) { /*Loop over rows to get the implicit scaling information. */
    big=0.0; 
    for (j=1;j<=n;j++) 
      if ((temp=fabs(a[i][j])) > big) big=temp;
    if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
    /*No nonzero largest element. */
    vv[i]=1.0/big; /*Save the scaling.*/
  }
  for (j=1;j<=n;j++) { /*This is the loop over columns of Crout's method.*/
    for (i=1;i<j;i++) { 
      sum=a[i][j]; 
      for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum; 
    } 
    big=0.0; /*Initialize for the search for largest pivot element.*/
    for (i=j;i<=n;i++) { 
      sum=a[i][j]; 
      for (k=1;k<j;k++)
	sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
      if ( (dum=vv[i]*fabs(sum)) >= big) { 
	/*Is the figure of merit for the pivot better than the best so far? */
	big=dum; 
	imax=i; 
      } 
    } 
    if (j != imax) { /*Do we need to interchange rows?*/
      for (k=1;k<=n;k++) { /*Yes, do so... */
	dum=a[imax][k];
	a[imax][k]=a[j][k];
	a[j][k]=dum; 
      } 
      *d = -(*d); /*...and change the parity of d.*/
      vv[imax]=vv[j]; /*Also interchange the scale factor.*/
    } 
    indx[j]=imax;
    if (a[j][j] == 0.0) a[j][j]=TINY; 
    /*If the pivot element is zero the matrix is singular (at least to the precision
      of the algorithm). For some applications on singular matrices, it is desirable 
      to substitute TINY for zero. */
    if (j != n) { 
      /*Now,  nally, divide by the pivot element. */
      dum=1.0/(a[j][j]);
      for (i=j+1;i<=n;i++) a[i][j] *= dum; 
    }
  } /*Go back for the next column in the reduction. */
  free_vector(vv,1,n); 
}

/* ------------------------------------------------------------ */

void sor(double **a, double **b, double **c, double **d, double **e, double **f, double **u, int jmax, double rjac) {

  /*Successive overrelaxation solution of equation (19.5.25) with Chebyshev acceleration.
    a, b, c, d, e, and f are input as the coefficients of the equation, each dimensioned
    to the grid size [1..jmax][1..jmax]. u is input as the initial guess to the solution,
    usually zero, and returns with the final value. rjac is input as the spectral radius
    of the Jacobi iteration, or an estimate of it. */

  void nrerror(char error_text[]); 
  int ipass,j,jsw,l,lsw,n; 
  double anorm,anormf=0.0,omega=1.0,resid; 
  /*Double precision is a good idea for jmax bigger than about 25.*/

  for (j=2;j<jmax;j++) 
    /*Compute initial norm of residual and terminate iteration when norm has been reduced by a factor EPS.*/
    for (l=2;l<jmax;l++) anormf += fabs(f[j][l]); /*Assumes initial u is zero.*/
  for (n=1;n<=MAXITS;n++) { 
    anorm=0.0; 
    jsw=1; 
    for (ipass=1;ipass<=2;ipass++) { /*Odd-even ordering.*/ 
      lsw=jsw; 
      for (j=2;j<jmax;j++) { 
	for (l=lsw+1;l<jmax;l+=2) { 
	  resid=a[j][l]*u[j+1][l] +b[j][l]*u[j-1][l] +c[j][l]*u[j][l+1] +d[j][l]*u[j][l-1] +e[j][l]*u[j][l] -f[j][l];
	  anorm += fabs(resid); 
	  u[j][l] -= omega*resid/e[j][l]; 
	} 
	lsw=3-lsw; 
      } 
      jsw=3-jsw; 
      omega=(n == 1 && ipass == 1 ? 1.0/(1.0-0.5*rjac*rjac) : 1.0/(1.0-0.25*rjac*rjac*omega)); 
    } 
    if (anorm < EPS*anormf) return; 
  } nrerror("MAXITS exceeded"); 

}

/* ------------------------------------------------------------------------- */


void gaussj(double **a, int n, double **b, int m) {
  /*Linear equation solution by Gauss-Jordan elimination, 
    equation (2.1.1) above. a[1..n][1..n] is the input matrix. 
    b[1..n][1..m] is input containing the m right-hand side vectors. 
    On output, a is replaced by its matrix inverse, and b is replaced 
    by the corresponding set of solution vectors. */

  int *indxc,*indxr,*ipiv; int i,icol,irow,j,k,l,ll; 
  double big,dum,pivinv; 

  indxc=ivector(1,n); 
  /*The integer arrays ipiv, indxr, and indxc are used for bookkeeping on the pivoting. */
  indxr=ivector(1,n); 
  ipiv=ivector(1,n); 
  for (j=1;j<=n;j++) ipiv[j]=0; 
  for (i=1;i<=n;i++) { /*This is the main loop over the columns to be reduced. */
    big=0.0; 
    for (j=1;j<=n;j++) 
      /*This is the outer loop of the search for a pivot element. */
      if (ipiv[j] != 1) 
	for (k=1;k<=n;k++) { 
	  if (ipiv[k] == 0) { 
	    if (fabs(a[j][k]) >= big) { 
	      big=fabs(a[j][k]); 
	      irow=j; 
	      icol=k; 
	    } 
	  } 
	  else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1"); 
	} 
    ++(ipiv[icol]); 
    /*We now have the pivot element, so we interchange rows, if needed, 
      to put the pivot element on the diagonal. The columns are not physically 
      interchanged, only relabeled: indxc[i], the column of the ith pivot element,
      is the ith column that is reduced, while indxr[i] is the row in which that 
      pivot element was originally located. If indxr[i] 6 = indxc[i] there is an 
      implied column interchange. With this form of bookkeeping, the solution b's 
      will end up in the correct order, and the inverse matrix will be scrambled 
      by columns. */
    if (irow != icol) { 
      for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l]); 
      for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l]); 
    } 
    indxr[i]=irow; 
    /*We are now ready to divide the pivot row by the pivot element, located at irow and icol. */
    indxc[i]=icol; 
    if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2"); 
    pivinv=1.0/a[icol][icol]; 
    a[icol][icol]=1.0; 
    for (l=1;l<=n;l++) a[icol][l] *= pivinv; 
    for (l=1;l<=m;l++) b[icol][l] *= pivinv;
    for (ll=1;ll<=n;ll++) /*Next, we reduce the rows...*/
      if (ll != icol) { /*...except for the pivot one, of course.*/
	dum=a[ll][icol]; 
	a[ll][icol]=0.0; 
	for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum; 
	for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum; 
      } 
  } 
  /*This is the end of the main loop over columns of the reduction. 
It only remains to unscram- ble the solution in view of the column interchanges. 
We do this by interchanging pairs of columns in the reverse order that the permutation 
was built up. */
  
  for (l=n;l>=1;l--) {
    if (indxr[l] != indxc[l]) 
      for (k=1;k<=n;k++) SWAP(a[k][indxr[l]],a[k][indxc[l]]); 
  } 
  /*And we are done.*/
  free_ivector(ipiv,1,n); 
  free_ivector(indxr,1,n); 
  free_ivector(indxc,1,n); 
}

/* -------------------------------------- */
/* -------------------------------------- */

double Bernoulli ( double x) {

  double simple = 0;
  double x1,x2,x3,x4,x5;
   
  /*  X1,X2,X3,X4,X5 are calculated with the program Selberherr.f */


  x1 = -1000.0;
  x2 = -7.450580596923828e-08;
  x3 = 7.450580596923828e-08;
  x4 = 38.0; 
  x5 = 1038.0;

  if(simple) {
    if(fabs(x)<1.0e-10) return 1.0;
    else return x/(exp(x)-1.0);
  }
  else {
    if(x <= x1) return -x;
    else if(x>x1 && x<x2) return x/(exp(x)-1.0);
    else if(x>x2 && x<=x3) return 1.0 - x/2.0;
    else if(x>x3 && x<x4) return (x*exp(-x))/(1.0-exp(-x));
    else if(x>=x4 && x<x5) return x*exp(-x);
    else  return 0.0;
  }

}

/* -------------------------------------- */
/* -------------------------------------- */

/*  In(x)Ga(1-x)N This routine evaluates the chemical composition
    of InGaN alloy when are given the two lattice constants as given
    for example from x-ray difraction. It takes into account the fact
    that there is strain. */
double chem_composition(double cl, double al) {

  double P, QQ, R ,S, left, right, x, dy, y;
  double ca,cb,aa,ab,va,vb;
  // FILE * fp;

  ca = 5.705; // InN c-lattice constant (in amstrong)
  cb = 5.187; // GaN c-lattice
  aa = 3.533; // InN a-lattice
  ab = 3.185; // GaN a-lattice
  va = 0.32; // InN Poisson ratio v=2C13/C33
  vb = 0.3; // GaN Poisson ratio

  P = (va-vb)*(aa-ab)*(ca-cb);
  QQ = (1.0+vb)*(aa-ab)*(ca-cb)+(va-vb)*((aa-ab)*cb+(ca-cb)*(ab-al));
  R = (aa-ab)*((1.0+vb)*cb-cl)+(ca-cb)*((1.0+vb)*ab-vb*al)+(va-vb)*(ab-al)*cb;
  S = (1.0+vb)*ab*cb-vb*al*cb-ab*cl;

  /*
  fp = fopen("OUTPUT_FILES/func.dat","w");
  x = -1;
  while(x<2) {
    fprintf(fp,"%le  %le\n",x,S + ( R + ( QQ + P*x)*x)*x);
    x += 0.01;
  }
  fclose(fp);
  */

  left = S;
  right = P+QQ+R+S;

  //  if((left>0 && right >0) || (left<0 && right <0)) nrerror("1 Chem_position: no root possible");

    x = 0.5;
  do {
    dy = R + (2.0*QQ+3.0*P*x)*x;
    y = S + ( R + ( QQ + P*x)*x)*x;
    x -= y/dy;
  } while (fabs(y/dy)>1e-4 && x>0 && x<1);

  if(x>=1 || x<=0) nrerror("2 Chem_position: no root possible");

  return x;
  

}


//****************************************************************

/*	This function calculates eigenvalues and eigenfunctions for a square well in the
	general case of an heterojunction ( massqw is different from mass0).
	V is the potential barrier.  xw is the width of the well and x0 is the width of
	the barrier. In *file1 the routine dumps the even eigenfunctions, in *file2 the odd ones.
	x0 and xw are expressed in nanometers, V in eV. The masses ar relative.
	ESTEP = energy steps
*/

void analytical(double V, double x0, double xw, double mass0, double massqw, char * file1, char * file2) {

  const int max	= 10;
  const double ESTEP = 100.0;

  int j,odd_count=0,even_count=0;
  double E,k1,k2,k1_next,k2_next,m_odd[max],m_even[max],h,coeff,d;
  double next_val_odd,next_val_even,new_val_odd,new_val_even,old_val_odd,old_val_even;
 
  FILE * data1, * data2;

  fopen_s(&data1,file1,"w");
  if (!data1) {
    printf("Error open %s.\n",file1);
    exit(1); 
  }
  fopen_s(&data2,file2,"w");
  if (!data2) {
    printf("Error open %s.\n",file2);
    exit(1);
  } 
 
  h=1.0/ESTEP;
  E=h;		
  k1_next = sqrt(mass0*KKK*(V-E));
  k2_next = sqrt(massqw*KKK*E);
  old_val_even = 100;
  old_val_odd = 100;
  do {
			
    k1 = k1_next;
    k2 = k2_next;
    k1_next = sqrt(mass0*KKK*(V-(E+h)));
    k2_next = sqrt(massqw*KKK*(E+h));

    new_val_even = fabs(-k2*tan(k2*xw/2.0)+k1);
    next_val_even = fabs(-k2_next*tan(k2_next*xw/2.0)+k1_next);
    if(new_val_even<=old_val_even && new_val_even<next_val_even) {
      m_even[even_count]=E;
      even_count++;
    }
    old_val_even=new_val_even;

    new_val_odd=fabs(k2*(1.0/tan(k2*xw/2.0))+k1);
    next_val_odd=fabs(k2_next*(1.0/tan(k2_next*xw/2))+k1_next);
    if(new_val_odd<old_val_odd	&& new_val_odd<next_val_odd) {
      m_odd[odd_count]=E;
      odd_count++;
    }
    old_val_odd=new_val_odd;

    E=E+h;
  } while(E<V);

  for(j=0;j<even_count;j++)	printf(" Even= %lf\n",m_even[j]);
  for(j=0;j<odd_count;j++)	printf(" Odd= %lf\n",m_odd[j]);

  d=-x0-xw/2;
  while(d<-xw/2) {
    fprintf(data1,"%le  ",d);
    for(j=0;j<even_count;j++)	{
      k1 = sqrt(mass0*KKK*(V-m_even[j]));
      k2 = sqrt(massqw*KKK*m_even[j]);
      coeff = cos(-k2*xw/2.0)*exp(k1*xw/2.0);
      fprintf(data1,"%le  ",coeff*exp(k1*d));
    } 
    fprintf(data1,"\n");
    d=d+1.0e-2;
  }
  while(d<xw/2.0) {
    fprintf(data1,"%le  ",d);
    for(j=0;j<even_count;j++) {
      k2 = sqrt(massqw*KKK*m_even[j]);
      fprintf(data1,"%le  ",cos(k2*d));
    }
    fprintf(data1,"\n");
    d=d+1.0e-2;	
  }
  while(d<x0+xw/2.0) {
    fprintf(data1,"%le  ",d);
    for(j=0;j<even_count;j++) {
      k1 = sqrt(mass0*KKK*(V-m_even[j]));
      k2 = sqrt(massqw*KKK*m_even[j]);
      coeff = cos(k2*xw/2)*exp(k1*xw/2);
      fprintf(data1,"%le  ",coeff*exp(-k1*d));
    }
    fprintf(data1,"\n");
    d=d+1.0e-2;	
  }

  d=-x0-xw/2;
  while(d<-xw/2) {
    fprintf(data2,"%le  ",d);
    for(j=0;j<odd_count;j++) {
      k1 = sqrt(mass0*KKK*(V-m_odd[j]));
      k2 = sqrt(massqw*KKK*m_odd[j]);
      coeff = sin(-k2*xw/2)*exp(k1*xw/2);
      fprintf(data2,"%le  ",coeff*exp(k1*d));
    }
    fprintf(data2,"\n");
    d=d+1.0e-2;	
  }
  while(d<xw/2.0) {
    fprintf(data2,"%le  ",d);
    for(j=0;j<odd_count;j++)	{
      k2 = sqrt(massqw*KKK*m_odd[j]);
      fprintf(data2,"%le  ",sin(k2*d));
    }
    fprintf(data2,"\n");
    d=d+1.0e-2;
  }
  while(d<x0+xw/2.0) {
    fprintf(data2,"%le  ",d);
    for(j=0;j<odd_count;j++) {
      k1 = sqrt(mass0*KKK*(V-m_odd[j]));
      k2 = sqrt(massqw*KKK*m_odd[j]);
      coeff = sin(k2*xw/2)*exp(k1*xw/2);
      fprintf(data2,"%le  ",coeff*exp(-k1*d));
    }
    fprintf(data2,"\n");
    d=d+1.0e-2;		
  }			

  fclose(data1);
  fclose(data2);
}

// -----------------------------------
// -----------------------------------
// type = 0; -> random;
// type = 1; -> sin()-shape;

double * fluctuation(int type, int num, double mean, double perc ) {

  int i,count;
  double deviation = (perc/100.0)*mean,sum;
  double * tmp = new double [num];
  FILE * fp;

  if(type==0) {
    count = 0;
    do {
      srand(11051973);
      sum = 0;
      for(i=0;i<num;i++) {
	tmp[i]=mean+deviation*2.0*(((double)rand()/RAND_MAX)-0.5);
	sum += tmp[i];
      }
      sum = sum/(double)num;
      count ++;
    } while(count<10 && fabs(sum-mean)>1e-2);
  }
  else if(type==1) {
    count = 0;
    do { 
      sum = 0;
      for(i=0;i<num;i++) {
	tmp[i]=mean+deviation*sin((double)i*PI/15.0);
	sum += tmp[i];
      }
      count ++;
    } while(count<10 && fabs(sum-mean)>1e-2);
  }
  else nrerror("FLUCTUATION: wrong type!");


  fopen_s(&fp,"OUTPUT_FILES/fluct.dat","w");
  fprintf(fp,"#Sum=%le  mean=%le  diff=%le\n",sum,mean,fabs(sum-mean));
  for(i=0;i<num;i++) fprintf(fp,"%d  %le\n",i,tmp[i]);
  fclose(fp);

  return tmp;

}
